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Abstract 

A first study of numerical Monte Carlo simulations with two quark doublets, a 
mass-degenerate one and a mass-split one, interpreted as u, d, s and c quarks, is 
carried out in the framework of the twisted mass Wilson lattice formulation. Tuning 
the bare parameters of this theory is explored on 12^ • 24 and 16^ • 32 lattices at lattice 
spacings a ~ 0.20 fm and a ~ 0.15 fm, respectively. 
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1 Introduction 



In QCD the effect of virtual quark loops is most important for the three light quarks 
{u,d,s). In recent unquenched numerical simulations, besides the two lightest quarks 
u and d, the s-quark is also included (see, for instance, [D El)- The formulation of 
QCD with twisted-mass Wilson fermions (Hj is based on chiral rotations of the bare 
mass (or, equivalently, of the Wilson-term) within quark doublets. Therefore, in this 
formulation there are two possibilities for unquenched simulations with (n, d, s) quarks: 
either the s quark is taken alone and the twisted mass formulation is restricted to the 
(n, d) doublet or, in addition to the s quark, also the c-quark is included in a mass-split 
doublet using the formulation in Ref. 4 . In the present paper we explore the latter 
possibility (for first results along this line see the proceedings contribution 5 ). 

Numerical simulations with twisted-mass Wilson fermions are usually performed at 
(or near) the critical (untwisted) bare quark mass, because there an automatic 0{a) 
improvement of the continuum limit is expected [H]. Our collaboration has performed 
several studies of twisted-mass QCD both in the quenched approximation [7j, [S], 05 
[in] and in unquenched simulations with dynamical {u,d) quarks ^21, ^Bji 
In the present paper we explore the possibility of numerical simulations of QCD with 
a degenerate doublet {u, d) and an mass-split doublet (c, s) of dynamical quarks in the 
twisted mass Wilson formulation. 

The plan of this paper is as follows: in the next section we define the lattice action 
and describe the simulation algorithm. Section |21 is devoted to the introduction of 
physical fields and currents important for the interpretation of results. In Section 0] we 
present our numerical results. The last section contains a discussion and final remarks. 



2 Lattice action and simulation algorithm 
2.1 Lattice action 

The notation for the lattice action of the light mass-degenerate {u, (i)-doublet, denoted 
by a subscript /, is the same as in our previous papers, for instance, Ref. |14j : 

-S*/ = XI I + ^75^-3 afii]xi,x) {^i,x+fiUxt,[r + 7p,]xi,x) | 

= T.Xi,xQExi,y . (1) 



,,xy 



Here, and in most cases below, colour-, spinor- and isospin indices are suppressed. For 
the isospin indices later on we shall also use notations as, for instance, xi = iXu Xd)- 
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The ( "untwisted" ) bare quark mass of the hght doublet in lattice units is denoted by 

fii^i = amoi + 4r = , (2) 
Zki 

r is the Wilson-parameter, set in our simulations to r = 1, amoi is another convention 
for the bare quark mass in lattice units and m is the conventional hopping parameter. 
The twisted mass of the light doublet in lattice units is denoted by afn. Ux^ G SU(3) 
is the gauge link variable and we also defined Ux-^ = f^i„^ ^ and = —7^. 

Besides the quark doublet fields xuXi in © it will turn out convenient to introduce 
other fields by the transformation 

IpUx = -^{^ + ilbT-i) XUx , V'/.x = X/,x-^(l+n5T3) • (3) 

The quark matrix on the x-basis Q;'^^ defined in (© is 

= ^xy ifJ-Kl + il5T3 ^x,y+f,Uyf,[r + 7^] (4) 

^ ti=±l 

or in a short notation, without the site indices, 

<5{^^ = f^Ki + «75T3 afii + N + R, (5) 

with 

1 ±4 ±4 
^xy = ~9 ^x,y+p.Uyfj,Jn , Rxy = ~q 5x,y+jjUy^ . (6) 

On the ip -basis defined in (jSJ we have the quark matrix 

= ^ (1 - i75T3) Qi""^ (1 - i75T3) =a^il + N- i-f5T3 {fl^l + R) . (7) 



As it has been shown by Frezzotti and Rossi in Ref. [1], a real quark determinant 
with a mass-split doublet can be obtained if the mass splitting is taken to be orthogonal 
in isospin space to the twist direction. One could take, for instance, the mass splitting 
in the ti direction if the twist is in the direction, as in (©. It is, however, more 
natural to exchange the two directions because then the bare mass is diagonal in 
isospin. In this case, the lattice action of the heavy mass-split (c, s)-doublet, denoted 
by a subscript /i, is defined as 

Sh = YXh,xQh;lyXh,y (8) 

x,y 

with 

Q^h^ = fJ-Kh + «75n aiJ.a + Tsaj^s -\-N + R . (9) 
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For the isospin indices later on we shall also use notations as, for instance, Xh = {Xc Xs)- 
The ip -basis is introduced similarly to © by 

iph,x = -^{'i- + h5Ti)xh,x, i^h,x = Xh,x^{'^ + ii5n) (10) 

and the quark matrix on the ip -basis is for the heavy mass-split doublet 

g}f ) = i (1 - ij5Ti) Q^^^ (1 - i75ri) = ail, + rga/i^ + N- i^^n {fi^h + R) ■ (H) 

For the SU(3) Yang-Mills gauge field we apply the tree-level improved Symanzik 
(tlSym) action which belongs to a one-parameter family of actions obtained by renor- 
malisation group considerations and in the Symanzik improvement scheme |15J . Those 
actions also include, besides the usual (1 x 1) Wilson loop plaquette term, planar 
rectangular (1x2) Wilson loops: 

^5 = /3E(co E |l-^ReC/i;j|+ci ii-\^R^ul;^\\ (12) 

with the condition cq = 1 — 8ci. For the tlSym action we have ci = —1/12 ll(ij. 



2.2 Simulation algorithm 

For preparing the sequences of gauge configurations a Polynomial Hybrid Monte Carlo 
(PHMC) updating algorithm was used. This algorithm is based on multi-step (actually 
two-step) polynomial approximations of the inverse fermion matrix with stochastic 
correction in the update chain as described in Ref. jJTj. It is based on the PHMC 
algorithm as introduced in Ref. jlSj . The polynomial approximation scheme and the 
stochastic correction in the update chain is taken over from the two-step multi-boson 
algorithm of Ref. (For an alternative updating algorithm in QCD with Nj = 

2 + 1 + 1 quark flavours, which will be used for algorithmic comparisons in the future, 
see 120] ■) 

For typical values of the approximation interval and polynomial orders on 16'^ • 32 
lattices see Table ^ The notations are those of Ref. : the approximation interval is 
[e. A], the orders of the polynomials Pj (j = 1, 2) are nj and those of Pj [j = 1, 2) are 
fij, respectively. The simulations have been done with determinant break-up ub = 2. 
On the 12^ • 24, for similar values of the pseudoscalar masses in lattice units, the orders 
n2 and n2 are the same and the values of ni and ni are somewhat smaller. 



3 Physical fields and currents 

The physical quark fields, which in the continuum limit are proportional to the renor- 
malised quark fields of both flavours in the doublets, are obtained 0| by a chiral 
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rotation from the fields in the lattice action in or from those defined in (jSJ for the 
light doublet, and similarly in ||H))- H1U() for the heavy doublet. On the x-basis we have 

^vhys ^ e^^^^-^^xi,. , '^S'' = Xi,.ei-'''^^ ; (13) 

<r = et-'^^-^XM, = XM^^"""^" ■ (14) 

Since the transformations in Q and pOj) correspond to chiral rotations with ui = ^ 
and LOh = respectively, we have with 

_ 71" _ vr 

f^i = - 2 ' = w/, - - (15) 

the relations 

V^,,,ei^'^«-3 . (16) 

- i^h,.e^'"'^^^' ■ (17) 

Since the simulations are usually performed near full twist corresponding to uji = uJh = 
^, the modified twist angles are close to zero: 

uJl^O, OhC^O . (18) 

Therefore, near full twist the ^-fields are approximately equal to the physical quark 
fields. At full twist the use of the V-basis is advantageous because the formulas are 
simpler than in the x-basis. 

The definition of the twist angles is not unique. There are different viable possibil- 
ities to define them and the critical hopping parameters corresponding to them (see, 
for instance, jH], [Hj, [HI, |1S|, PH, US). 

Here, for the light doublet, we use the definition based on the requirement of parity 
conservation for some matrix element of the physical vector and axialvector current, 
as first introduced in |2^, jT2] and numerically studied in detail in 14 . For this let us 
introduce the bare vector- and axialvector bilinears 

Kx), = Xi,x\ra7^xi,x , Al^^ = Xi,x\Taliilr>XUx {a = 1,2). (19) 

The twist angle is introduced as the chiral rotation angle between the renormalised 
(physical) chiral currents: 

yl%^. = ZivVi^^cosui +eabZiAAl,^s\nuji , (20) 
Alx^, = ZiAAf.,^ cos uJi +eabZivVi'',.^ sin uJi (21) 

where only charged currents are considered (a=l,2), Cab is the antisymmetric unit 
tensor and Ziy and Zi^ are the multiplicative renormalisation factors of the vector 
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and axialvector current, respectively. The exact requirements defining lvi (and also 
yielding the value of Zia/Ziy) is taken to be 



( I ^=0 ) = , ( I it,,^=i,2,3 I ) = . (22) 

For the heavy doublet, in principle, one could translate and use this construction, 
too, but for applications in the kaon and D-meson sector it is more natural to consider 
bilinears between the light and the heavy doublet. In addition, inside the heavy dou- 
blet, due to the off-diagonal twist, one also would have to consider disconnected quark 
contributions which are absent in the light-heavy sector. Let us introduce the bare 
vector-, axialvector-, scalar- and pseudoscalar bilinears in the K~^- and D'^-sector as 

(23) 
(24) 
(25) 
(26) 

and similarly for the K^- and L>~-sector by changing u ^ d. Denoting the kaon- and 
D-meson-doublet hy K = {K^ K^) and D = {D^ D~), respectively, and introducing 
K = {K~^ —K^) and D = [D^ —D^), the renormalised vector and axialvector currents 
of the kaon doublet are given by 

VK,x^l = cos -- cos —ZvVK,x^l + sm -- sm —ZyV^ 



^K+,xfi 


~ Xs,x1 iJ.Xu,x ) 


^K+,xii 


= Xs,xlf^75Xu,x 


Sk+,x 


= Xs,xXu,x 1 


Pk+,x 


= Xs,x75Xu,x , 


^DO,xii 


= Xc,x^fJ.Xu,x ) 


^DO,xfi 


= Xc,xltilbXu,x 


Sdo,x 


= Xc,xXu,x ) 


Pdo,x ^ 


= Xc,x75Xu,x , 



— cos — ZAAi^. —I cos — sm — 
2 2 ^'^'^ 2 2 

D,xii 



+ tsm — cos —ZaA^^^^ - z cos y sm —ZaAd,x/i , (27) 
AK,x^l = cos — cos —ZaAk,xii + sm — sm —ZaA 



I sm — cos —ZvVf. ^ - z cos — sm — Zy Vd.x/^ ■ (28) 



2 2 ^'^^ 2 2 

Analogously for the scalar bilinears: 



^K,x = cos — cos —Zst>K,x - sm — sm —ZsS^^^ 



z sm y COS —ZpPp.^ + ICOS — sm —ZpPd,x , (29) 



Pk,x = COS y COS —ZpPk,x " sm — sm —ZpP^ ^ 



+ I sui — COS— Zsb p. + I COS — sm.—ZsbD,x ■ (30) 

Similar relations hold in the D-meson doublet, too. (|27|) - (|28j) show that near full twist 
^l,h — vr/2 all four terms on the right hand sides have roughly equal coefficients. 
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The requirement of parity symmetry in the isotriplets (pions, rho-mesons) ahows to 
fix the twist angle loi, cf. H22|) . In the case of the heavy-hght isodoublet one has to take 
into account the mixing between the kaons and D-mesons. In this case the twist angle 
LVh (and uJi) can be fixed by requiring conservation of parity and/or flavour symmetry. 

The equations in (|22j) follow by considering |2j, |12j . |14j the (vanishing) vector- 
current-pseudoscalar and axialvector-current- vector-current correlators, which turns 
out to be the most convenient choice for fixing the twist angle in the light sector. In 
the heavy-light sector the mixing patterns for currents and scalar bilinears are similar, 
so any combination of operators gives similar expressions. However, correlators only 
made up of scalar bilinears are expected to give a better signal, so we concentrate on 
this case for the discussion. Considering the upper components, four bilinears Pk+, 
Pf)0, Sx+, Sdo and the respective charge-conjugated versions must be included in 
the analysis. We define a four-dimensional vector of the multiplicatively renormalised 
bilinears 



V 



f ZpPk+\ 
ZpPoo 
ZsSk+ 

\ZsSdo ) 



V = {-ZpPx-,-ZpP£)0,ZsSK-,ZsS 



DO J 



(31) 



and analogously the vector V of the fully renormalised bilinears according to ()29|). 
(|3()|) (and the analogous equations for the Z)-mesons). and (jHU]) can be then 

reformulated in a compact notation as 



V = MV 



V = VM' 



(32) 



with the 4x4 matrix A4 given by 
/ 



M{uJi,uJh) 



V 



T 


cos^f 


-sin^f 


sin^f 


i sin ^ 


cos^ 


isin^^ 


cosf 


-2 


sin^f 


cos^f 


cos^ 


isin^f- 


cos^ 


i sin ^ 


cos^^ 


2 


cos^^ 


i sin ^ 


cosf 


cos^f 


cos^^ 


-sin^f 


sin^f 


2 


cos^f 


i sin ^ 


cos^ 


-sinf 


sin^^ 


cos^f 


cos^^ 



(33) 

Ai is the unitary matrix describing the mixing pattern between the kaon and D-meson 
doublets. One can easily see that 

M'^ = M and M^LOi, ujh) = M* {uji, LOh) = M{-uJi, -LOh) = M-^ {ioi, iOh) . (34) 

(The last equality is expected since reversing the sign of the angles corresponds to the 
inverse chiral transformation). One can at this point define a correlator matrix in the 
kaon-D-meson sector by 

C = (V0V) (35) 



(for example, Cu = —Zp{Px+PK-)) and its fully renormalised version C = {V V). 
One has 



C = M{uJi,ujh)CM ^{uji,ujh) , 



(36) 
(37) 



Restoration of parity- and flavour-symmetry implies that C is a diagonal matrix with 
Ai{uji,uJh) the matrix realizing the diagonalisation. The off-diagonal elements of the 
matrix equation (|36|) can be in principle used to determine the angles uji and ujh, while 
the diagonal elements give the physical correlators from which e.g. the masses can be 
obtained. Of course, in general, parity and flavour can only be restored up to 0{a) 
violations. 

Taking also into account the residual discrete symmetries possessed by the action 
defined by ^ and ((SJ-®, the only non-trivial conditions are obtained by imposing the 
vanishing of the flavour violating matrix elements C12, C34 and transposed. Defining 



for brevity si = sin Sh = sin q = cos Ch = cos the conditions are: 



C12 + c 



21 



[ciCh)'^ + {siShf (Cl2+C2l)+ (siCh)'^ + {ShClf (C34+C, 



43 J 



-2ciChSiSh{Cii + C22 - C33 - C44) + ishCh{sf - cf){Ci3 



C31+C 



24 



42 J 



+iSlCl{sh^ - Ch'){C; 



23 



C32+C 



14 



--41 J 



, 



(38) 



C34 + C. 



43 



{ChSlf + {ciShf (Cl2+C2l)+ {ciChf + {siShf (C34+C, 



43 J 



+2ciChSiSh{Cii + C22 - C33 - C44) - ishCh{sf - )(Ci3 - C31 + C24 - C. 



42 J 



-iSlCl{sh^ - Ch^){C23 - C32 + Ci4 - C41) = . 

The sum of the two above equations implies 

C12 + C21 -I- C34 -I- C43 = . 



(39) 



(40) 



A non trivial relation for the renormalisation constants of the bilinears is obtained 
from (001) 



(41) 



{Pk+Pdo) + {PdoPk-) ' 
which can be used for a non-perturbative determination of Zp/Zs- 

Using (|in|) . (|nH|) (or (jBHI)) can be restated in a compact way as a relation between 
cot ujh and cot uJi 

Cu + C22 - C33 - C44 + i{Ci3 - C31 + C24 - C42) cot cvi 



cot UJh 



(C12 + C21 - C34 - C43) cot UJl - i{C23 - C32 + Ci4 - C. 



(42) 



41 



This can be used to determine u)h once u)i is known, {loi can be obtained following the 
prescription of [H], [H], HU). 



8 



This discussion suggests that, especiahy near full twist where the mixing is max- 
imal, the analysis of the masses in the kaon-D-meson sector should be performed by 
considering the 4-dimensional correlator matrix C. 

For tuning the hopping parameters the untwisted PC AC quark mass is also very use- 
ful. In the light doublet it is defined by the PCAC-relation containing the axialvector 
current Af^^ in ((TO)) and the corresponding pseudoscalar density = Xi,x'^'^a'y5Xi,x- 



If)* 4+ P 

^^PCAC ^ (43) 



where r± = ri it ir2. The condition of full twist in the light quark sector obtained 
from (|^^ by setting uoi = t: /2 coincides ^1] with m^f^'-^ = 0. 

In the heavy sector one can define an untwisted PCAC quark mass Jtz^^'^'^, too. 
A natural definition is obtained by considering the axialvector Ward identity 

2iZ]^^a^jLaSl,j,, a = l 
5*A^,,^ = 2am^f^^P,% + <j 0, ' a = 2 (44) 

(-2i)Z/a/i5P0,, a = 3 

where, in analogy with the light sector in (|19)) . we define 

'^h,xii = Xh,X2'^aliil5Xh,x {a = 1, 2, 3) , 5° 3, = Xh,xXh,x , Ph,x = Xh,xl5Xh,x ■ 

(45) 

(Observe that for uniformity with the definition l|43j) we incorporate a factor Z^^ in the 
definition of the untwisted PCAC quark mass). The above identity could in principle 
be used to tune ujh to -it/2 by imposing am^^^'-^ = 0. However, as already mentioned, 
the presence of disconnected contributions in the heavy sector are likely not to allow 
for precise determinations. 

One can consider also in this case the heavy- light sector. Here the axialvector Ward 
identities read 

d*AK,x,, = {am^^^^ + am^f'^^)PK,x^^ + iZJ^a^llSj^^^^ + iZJ^a^I„SD,x^. (46) 
5*^D,xM = {am^c^^ + am^f^^^) PD,xf, + iZj^ aw S^^^^ + iZ^^ a/x^ SK,x,^ ■ (47) 

The solution of the over-determined linear system, obtained by taking a suitable matrix 
element (for instance, {d*Aj^+ .j.^Pj^- y)), allows to determine numerically (together 
with (|43|) ) the untwisted PCAC mass of the heavy quarks m^^^*^, m^^^'^ and the 
renormalisation factor Z^. The condition of full twist in the heavy doublet can be 
written as 

mP^^C ^ ^PCAC ^ ^PCAC ^ _ ^4g) 

The quark masses defined by (|43|) and (|46j) - (|47j) are untwisted components of bare 
quark masses. The physical quark masses can be obtained by the corresponding PCAC- 
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relations of the renormalised currents and densities: 



PCAC I PCAC — i^tJ,^K+,xfi Pr- ,y) ,(.„n 

anig + amf = — ^ — , (50) 

{Pk+,x PK-,y) 
PCAC I PCAC — ^^fJ,^D+ ,xfjLPD- ,y) 

am^ + amf = — ^ h: — . (51) 

{Pd+,x PD-,y) 



They are related to the bare quark masses by 



m 



PCAC 

I - ^p yy^Ai'^i 



Zp'J{ZAmPf•^cr + ^,f, (52) 



<f'"' = Zp'^{ZAmP<^^cy + ^,l±Z,'^^s. (53) 

4 Numerical simulations 

Our main goal in this work is to gain experience with the tuning of lattice parameters 
for future large scale simulations. Based on our recent work with A'^^ = 2 dynamical 
twisted mass Wilson fermion QCD simulations in Refs. jJJ, JS]) !14j and j26j . 
the main emphasis is on the effects of the additional dynamical flavours s and c. As 
in the Nf = 2 case, we start with coarse lattices: lattice spacings about a ~ 0.2 fm 
on a 12^ • 24 lattice and a ~ 0.15 fm on a 16^ • 32 lattice. (This implies spatial lattice 
extensions of L ~ 2.4 fm.) The parameters of our main runs are: on the 12^ • 24 
lattice f3 = 3.25, a/x; = 0.01, a/Xo- = 0.315, afis = 0.285 and on the 16^ • 32 lattice 
(3 = 3.35, a/i/ = 0.0075, a^a = 0.2363, afis = 0.2138. The statistics is between 500 
and 1100 PHMC trajectories of length 0.4. (Of course, in order to find the appropriate 
parameters, we also had to perform at the beginning several additional short runs 
which we do not include here.) 

The tuning to full twist of the theory with an additional heavy doublet is com- 
plicated by the fact that two independent parameters m and Kh must be set to their 
respective critical values, using e.g. for the heavy sector the procedure outlined in the 
previous section. However, it can be shown that in the continuum limit the deviation of 
the two critical hopping parameters Ki^cr and Kh,cr goes to zero as 0{a). An argument 
is given in the Appendix. This suggests to tune ki to the value where m^f^^^ = with 
Kh = Ki- in this situation m^^^"-^ = 0{a). Observe that since the average quark mass 
in the heavy sector is typically large, the 0{a) error is expected not to affect the full 
twist improvement in the sense of ,6;, while it is critical to have good tuning in the 
light quark sector. This can be checked by computing ujh as suggested in the previous 
section and verifying uj^ ~ tt/2. 
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In view of this, we have set the two hopping parameters to be equal in our main 
runs: k = ki = hh- (In ^ few additional runs we checked that small individual changes 
of by Ak/i ~ 0.001 do not alter any of the qualitative conclusions.) 

The average plaquette values as a function of the hopping parameter m = Kh = k, 
for fixed values of the twisted masses, are shown by Figures lU and El on the 12^ • 24 
and 16^ • 32 lattices, respectively. On the 12^ • 24 lattice a strong metastability is 
observed for 0.1745 < k < 0.1747, which we interpret as the manifestation of a first 
order phase transition. This behaviour agrees with one of the scenarios predicted by 
ChPT including leading lattice artifacts [231, |2HI, [22], [201 • It has also been observed 
in our previous simulations, for instance, in JJ. On the 16^-32 lattice no metastability 
could be observed, although there is a sharp rise of the average plaquette value between 
n = 0.1705 and k = 0.1706. This may also signal a (weaker) first order phase transition 
or a cross-over. To decide among these two possibilities, in principle, an investigation 
of the infinite volume behaviour would be necessary. In practice, in a finite volume, 
the effects of a real first order phase transition and a cross-over are similar. 

We emphasize that this observed behaviour is not related to some imperfection of 
the simulation algorithm. Due to the positive twisted masses the eigenvalues of the 
fermion matrix have a positive lower bound. Therefore, we could choose the HMC step 
size small enough in order that the molecular dynamical force does not become too 
large. The behaviour of the system when crossing the phase transition region is nicely 
illustrated by the run history in Figure El One can recognize three stages in the plot: 
metastable start at m^f"^'^ > 0; crossing; stable thermalization at m^f^'-' < 0. A high 
concentration of small eigenvalues occurs during the crossing, because a large portion 
of the Dirac spectrum (actually all the physically relevant eigenvalues) is moving from 
the right half complex plane with Re A > to the left one with Re A < 0. 

We determined several quantities in both the pion- and kaon-sector. The values of 
some of them are collected in Tables |21 and 01 As in our previous work, we determined 
the lattice spacing from the quark force by the Sommer scale parameter assuming 
ro = 0.5 fm. Taking the values for positive untwisted PCAC quark masses {am^f^^'-"' > 
0), we get for /3 = 3.25 on the 12^ • 24 lattice a{(3 = 3.25) ~ 0.20 fm and for p = 3.35 
on the 16^ • 32 lattice a(/3 = 3.35) ~ 0.15 fm. These correspond to a^^ ~ 1.0 GeV and 
a^^ ~ 1.3 GeV, respectively. 

It follows from the data in Tables |2l and (HI that the pion, and hence the u- and 
d-quark masses, are not particularly small in our runs. Considering only the points 
with positive untwisted PCAC quark mass {am^f^'-" > 0) outside the metastability 
region at /? = 3.25 we have m-,^ > 670 MeV. At /? = 3.35 the corresponding value 
is > 450 MeV. (The points with am^p"^*^ < have > 530 MeV and > 
560 MeV for the two /3-values, respectively, but they are usually not considered for 
large scale simulations because of the strongly fluctuating small eigenvalues as shown. 
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for instance, by Fig. 13) 

The kaon masses are also given Tables HI and |21 Let us note that, in the Frezzotti- 
Rossi formulation of the split-mass doublet we use, the masses in the kaon doublet 
(and D-meson doublet) are exactly degenerate. This follows from an exact symmetry 
of the lattice action defined in Section 12.11 (both in the x- ^-iid '(/'-basis of quark fields) 
namely, simultaneous multiplication by an isospin matrix and space reflection: 



This exact symmetry exchanges the u-quark and the d-quark, hence the equality of the 
masses within kaon- and D-meson doublets follows. 

Let us note that in a recent publication |32j a non-zero kaon mass splitting has 
been calculated in the quenched approximation using another formulation [33j of the 
mass-split doublet where both the twist and the mass splitting are in the same isospin 
direction. This formulation has, however, the disadvantage that the fermion determi- 
nant is non-real and therefore an unquenched computation is practically impossible 
at present. The difference in the presence and absence of the kaon mass splitting in 
the two formulations comes from the fact that the states with a given quark flavour 
correspond to different linear combinations here and there. 

Similarly to the pion masses, the kaon masses in Tables [21 and 01 are also higher than 
the physical value. In the points cited above for the pion mass we have: at /3 = 3.25 
and 13 = 3.35 rriK > 920 MeV and rriK > 850 MeV, respectively. The kaon mass can be 
easily lowered by tuning the mass parameters in the heavy doublet. In order to explore 
this we also performed simulations at /? = 3.35, am = 0.0075 on the 16^ • 32 lattice 
with a/Xo- = afis = 0.15. For instance, at ki = = 0.17 we got am,r = 0.4432(40) 
and arriK = 0.5918(22). Comparing to the third line in Table El one can see that both 
the pion and the kaon mass become smaller. In particular, the kaon mass is smaller 
by a factor of about 3/4. This shows that the kaon mass can probably be tuned to 
its physical value if wanted. Another possibility is to do the chiral extrapolation by 
fixing, instead of rriK, the pion-kaon mass ratio m^/mx to its physical value. 

The D-meson masses in Tables El and O are typically smaller than the physical value. 
In the points cited above for the pion and kaon masses we have: at /? = 3.25 and 
/3 = 3.35 mo ~ 1450 MeV and rriD ~ 1400 MeV, respectively, can, in principle, 
also be tuned to its physical value. However, on coarse lattices the D-meson mass 
is close to the cut-off and, therefore, it is more reasonable to keep it smaller than 
the physical value in order to be well below the cut-off. In fact, in our runs the 
actual D-meson masses are already at the cut-off because we have ~ 1 GeV and 




(54) 
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a^^ ~ 1.3 GeV at /3 = 3.25 and /5 = 3.35, respectively. But on a fine lattice, say with 
Q-i ~ 4 GeV, it will become possible to directly go to the physical value of rriD, too. 

The machinery for the twist angle in the heavy doublet uj^ developed in Sec. 01 
has been tested in a few runs, too. The formulas worked fine and the results turned 
out to be plausible. For instance, in the run at /? = 3.35, ki = = 0.1704, a/i; = 
0.0075, a/Xg- = 0.2363, a/x^ = 0.2138 on a 16^ • 32 lattice we obtained from 400 gauge 
configurations: 

Wi/TT = 0.0981(55), tj/,/7r = 0.490(25), 

Zp/Zcj = 0.5739(65), = 0.897(11), = 0.5490(12). (55) 

As one sees, ujh is rather close to -7r/2 even if oji is still far from it. This is a consequence 
of fia- ^ fJ-i- Using the relation (valid in the continuum) cot(tj/i)/ cot(LL';) = fii/fia and 
the value of uji given above, one would get w/j/vr = 0.468. The situation is very similar 
in the runs on a 12^ • 24 lattice, too. For instance, in the run with largest untwisted 
mass of Table El at ki = Kh = 0.1740/, we obtained: 

uji/n = 0.04298(34), w/^/vr = 0.4356(83), Zp/Zs = 0.581(11). (56) 

These results imply that putting the untwisted quark mass equal in the two sectors 
gives an elegant solution for tuning to full twist: one can just do the same as in the 
Nf = 2 case. Due to the large twisted component in the heavy sector, the tuning of 
w/i to 7r/2 is no problem at all: already at moderate values of uji, ujh is almost equal to 
7r/2. 

Let us finally mention that using Chiral Perturbation Theory (ChPT) formulas one 
can also extrapolate from our simulation points to smaller pion- and kaon-masses. As 
a simple example, let us take the squared pion-kaon mass ratio in lowest order (LO) 
ChPT: 

{m-^/niKf = ^ • (57) 

mud + rus 

In terms of our parameters we can set 

m„, = y^(Z^m^P^^)2 + ^2 ^ = ^{Za m^f ^^)2 + (^,)2 - ^ (58) 

where Zp/Zs is a fitted relative renormalisation factor. In our fits we set, for simplicity, 
Za = 0.897 from and we also assumed m^^^"-^ = m^f"'^'^, which corresponds to 
the assumption K^^cr = i^i,cr- The results for both (3 values are shown in FigurejHl Note 
that although these fits look rather good, clearly, the validity of Chiral Perturbation 
Theory in general has to be checked in further simulations at small values of a and m,r- 
It turns out that the fitted values of Zp/Zs are well below 1, namely Zp/Zs — 0.45, 
which implies that, as also directly shown by our simulation data, the kaon mass reacts 
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relatively weakly to the change of the bare quark mass difference parameter a/x^. The 
deviation of Zp/Zs obtained in the LO-ChPT fit from the values in (jSSj) and (|56() 
might be due to lattice artifacts or/and to the fact that in ()55 |) -(|56 |) no extrapolation 
to zero quark masses is performed. 

Note that the obtained values of Zp/Zs do not satisfy the bound derived in 0] 
which would ensure the positivity of the quark determinant, because in case of the 
(c, s)-doublet this bound is Zp/Zs > ("^c — rns)/{mc + rris) ~ 0.85. This means that 
there might be some gauge configuration where the determinant of the (c, s)-doublet 
is negative. However, such configurations have a very low probability and hence they 
practically never occur in Monte Carlo simulations. This is shown by the eigenvalues 
of the fermion matrix which never come close to zero: for the (c, s)-doublet in our 
simulations they always satisfy Xmin.h > 0.01. (This has to be compared to the minimal 
eigenvalues in the (u, (i)-doublet which only satisfy Xmin,i > 0.0001.) 

It is remarkable that the minimum value of the interpolated curves in Figure 1^1 are 
not far away from the physical value (mT^/niK)'^ — 0.082. This raises the interesting 
question, whether it would be possible to perform uncoventional chiral extrapolations 
from simulation data at fixed twisted masses. 

5 Discussion 

The main conclusion of the present paper is that numerical simulations of QCD with 
unquenched u, d, s and c quarks are possible in the twisted-mass Wilson formulation. 

The PHMC updating algorithm with multi-step polynomial approximations and 
stochastic correction during the update turned out to be effective even in difficult 
situations near a first order phase transition (or cross-over). The autocorrelations of 
the quantities given in Tables [21 and |21 are typically of 0(1) in number of PHMC- 
trajectories (most of the time of length 0.4), therefore, it is worth to analyse the gauge 
configurations after every trajectory. 

At /3 = 3.25 (lattice spacing a ~ 0.20 fm) on our 12^ • 24 lattice we observed strong 
metastabilities suggesting a first order phase transition. This agrees with one of the 
scenarios predicted by ChPT including leading lattice artifacts j2Zl , |2H| , 120] , [^01 and 
has been observed previously in several QCD simulations with Wilson fermions |34j . 
[SSI, ISni, CU- At /? = 3.35 (lattice spacing a ~ 0.15 fm) on our 16^ • 32 lattice the 
phase transition becomes weaker but is still visible as a strong cross-over region with 
fast changes in several quantities. Compared to Nj = 2 runs at similar lattice spacings 
the first order phase transition becomes stronger for Nf = 2 + 1 + 1. (This agrees with 
the early observations in |35j.) 

The smallest simulated pion mass in a stable point with positive untwisted PCAC 
quark mass {am^f^^'^ > 0) at /? = 3.25 (a ~ 0.20 fm) and /3 = 3.35 (a ~ 0.15 fm) is 
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TTiTT — 670 MeV and — 450 MeV, respectively. Our expectation based on the ChPT 
formulas and on our previous experience is that, for instance, on a • 48 lattice with 
a ~ 0.10 fm the minimal pion mass at a/i = 0.005 will be somewhere in the range 
270 MeV < < 300 MeV. This is because at vanishing twisted masses m™" is 

going to zero as 0{a) and for positive twisted mass the decrease is somewhat faster. 
(The lower value of the estimate corresponds to the minimum of the extrapolated curve 
in Figure ini) 

The kaon mass in the present simulations is higher than the physical value but 
can probably be properly tuned by changing the twisted mass parameters in the c- 
s doublet. The D-meson mass is smaller than the physical value (i.e. the c-s mass 
splitting is smaller than in nature) but this is reasonable on coarse lattices in order to 
stay with it below the cut-off. On finer lattices (say, with a ~ 0.05 fm) one can try 
to tune also the D-meson mass to its physical value. A possible difficulty in properly 
tuning the mass splittings in the c-s doublet can be caused by the relative insensitivity 
of the masses to the bare mass-splitting parameter afis- This may imply the necessity 
of some extrapolations in the mass ratios. 

In case of the c-s-doublet the mass splitting is rather large because the renormalised 
quark masses satisfy (rric — ms)/{mc + mg) ~ 0.85. Therefore it is important to take 
into account the mass splitting. For the u-d-doublet, well above the scale of u and d 
quark masses, the mass degeneracy can be considered as a good approximation, but 
even in this case we have in nature (m^ — r?T,u)/(mrf + rriu) ~ 0.28. Hence also there, 
on a long run, the problem of the quark mass splitting within the doublet has to be 
tackled. 

In summary, our experience in this paper is rather positive both for the twisted- 
mass Wilson fermion formulation and for the PHMC algorithm we are using. This 
opens the road for future large scale QCD simulations with dynamical u, d, s and c 
quarks. 
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Appendix 



In the Nf = 2 theory, one possible definition of the critical quark mass mocr(50) A*) is 



given by the vanishing of the PCAC quark mass Due to chirality breaking 



the latter gets shifted: 



m: 



PCAC 

X 



mo -a ^f{go,amo,afi) , (59) 



with / a dimensionless function. On the basis of the symmetry of the action under 
parity x (/i ^ —fi) one can show that the additive renormalisation of the quark mass 
is even in and analyticity in turn implies 



f{go, amo, a^i) = f{go, amo) + 0{fi a ) , (60) 

where f {go, amo) is the shift for ordinary Nf = 2 QCD without twisted mass term. So 
the twisted mass term in the action only produces an 0{a) effect on the quark mass 
(with go and mo held fixed): 



m 



PCAC 
X 



rriQ-a ^ f{go,amo) + 0{a) . (61) 



The above argument can be easily generalized to the Nf = 2 + 1 + 1 theory. Here 
one has to make a distinction between the two sectors: 



mi^f- = moi-a^fi{go,amoi,amoh,ani,afi„,ans), (62) 



PCAC 

xi 

xh^'~^ = rnoh-a~'^fh{9o,amoh,amoi,afi^,a^i,ai2s). (63) 



The functions /; and fh are in this case even in ^u/, fih and fj,s^: similarly to A'^^ = 2, 
the associated terms in the action only affect the additive renormalisation of the quark 
mass by 0{a) terms. So we write: 



m^l 
m^f^ 



PCAC 



= moi-a ^f{go,amoi,amoh) + 0{a), (64) 
^•^^•^ - moh - a'^figo, amoh, amoi) + 0{a) , (65) 



where on the r.h.s. we have now the mass-shifts for the theory without twist and mass- 
splitting {Nf = 2 + 2 QCD): here the distinction between the two sectors is immaterial. 
From eqs. (|64|). (|65|) it follows immediately 



moi = moh = mo => m^f = m^f^^ + 0{a) . (66) 



^An additional symmetry in the heavy sector is needed for the argument, namely Xh,x exp {i'^TiJxh,. 
Xh.x Xh.,xeTq){-i§Ti} composed with ^5 -> -^5. 
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Tables 



Table 1: Algorithmic parameters in two runs on a 16'^ ■ 32 lattice at j3 = 
3.35, Ki = Kh = K,, afii = 0.0075, a/Xg- = 0.2363, afj,s = 0.2138 and with 
determinant break-up ub — 2. The first line for a given k shows the pion 
mass and the parameters for the light doublet, the second line the kaon 
mass and the parameters for the heavy doublet. 
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Tabic 2: Selected results of the runs on a 12^-24 lattice at f3 = 3.25, afii = 
0.01, aficr = 0.315, ans = 0.285. The subscript on k = ki = Kh denotes: 
L for "low" and H for "high" plaquette phase, respectively. 
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Table 3: Selected results of the runs on a 16^-32 lattice at ^ — 3.35, a//; = 
0.0075, ai^„ = 0.2363, a/^s = 0.2138. 
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Figure 1: The average plaquette on 12^ • 24 lattice at (5 — 3.25, a^i 
0.01, a/Xp. = 0.315, a^s — 0.285 as a function of k = ki — Kh- 
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Figure 2: The average plaquette on 16^ ■ 32 lattice at (3 — 3.35, a^i — 
0.0075, a/^o- = 0.2363, a^s — 0.2138 as a function of k = ki — Kh- 
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Figure 3: Run history on a 16^-32 lattice at 13 = 3.35, aiii = 0.0075, a/ia = 
0.2363, a/is — 0.2138, ki = Kh — 0.1706. This run started from a pre- 
vious one at K = 0.1705. On the horizontal axis the number of PHMC- 
trajectories (of length At — OA) is given. The average plaquette (upper 
curve, left scale ) and the smallest eigenvalue of the squared preconditioned 
fermion matrix Xmin (lower curve, right scale) are shown. The horizon- 
tal lines indicate the average plaquette after equilibration and the absolute 
minimum of Xmin, respectively. 
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Figure 4: The mass-squared ofpion and kaon as a function of the untwisted 
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Figure 5: The mass-squared ofpion and kaon as a function of the untwisted 
PCAC quark mass on a 16^ • 32 lattice at (5 — 3.35, a^i — 0.0075, a/io- = 
0.2363, = 0.2138. 
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Figure 6: Lowest order ChPT fit of the squared pion to kaon mass ratio 
as a function of the untwisted PC AC quark mass. Squares and triangles 
are data at (3 = 3.35, fii = 0.0075 for fi„ = 0.2363, fis = 0.2138 and 
/Xo- = 0.15, fis = 0.15, respectively. The fit to the formulas ([371)- ( T^) 
gives Zp/Zs = 0.446. Circles are data at (3 = 3.25, ni = 0.01, ^„ = 
0.315, = 0.285. The fit gives in this case Zp/Zs = 0.457. 
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